Simulación de variables aleatorias normales¶

Método de Box-Muller¶

Recordemos que si $X_1,\ldots, X_n \stackrel{\text{i.i.d}}{\sim} \text{Normal}(0,1)$ entonces $$\sum_{i=1}^n X_i^2 = Q \sim \chi^2_n,$$ donde $\chi^2_n$ denota una distribución Ji-cuadrada con $n$ grados de libertad.

Por otro lado $Y\sim \text{Exponencial}\left( \frac{1}{2} \right)$ coincide con la distribución $\chi^2_2$ por lo que $R^2=X_1^2 + X_2^2$ puede ser simulada con el método de la tranformación inversa para variables aleatorias exponenciales. Específicamente, se genera $R^2=-\log(U)/\left(1/2\right)=-2\log(U)$ haciendo uso de $U\sim \text{Uniforme}(0,1)$. La distribución de $X_1,X_2|R=r$ es uniforme en el círculo $\left\{ (x,y)\, : \, x^2 + y^2 =r \right\}$. Se sigue que

$$ Z_1 = \sqrt{-2\log(U_1)}\sin(2\pi U_2), \quad Z_2 = \sqrt{-2\log(U_1)}\cos(2\pi U_2) $$

cumple que $Z_1,Z_2\stackrel{\text{i.i.d}}{\sim} \text{Normal}(0,1)$.

Método de Marsaglia¶

Observamos que si $(V_1,V_2)$ es uniforme en $\left\{ (v_1,v_2)\, : \, v_1^2 + v_2^2 \leq 1 \right\}$ entonces $W=R^2=V_1^2+V_2^2 \sim \text{Uniforme}(0,1)$, $\Theta=\arctan\left(V_2/V_1 \right)\sim \text{Uniforme}(0,2\pi)$ y estas variables aleatorias son independientes.

Se sigue que $-2\log(W)\sim \text{Exponencial}(1/2)$ y $\left(\cos(\Theta),\, \sin(\Theta)\right)=\left( V_1/\sqrt{W}, \, V_2/\sqrt{W} \right)$ es uniforme en el círculo. Por el método de Box-Muller se tiene que

$$ Z_1 = \sqrt{-2\log(W)/W}V_1, \quad Z_2 = \sqrt{-2\log(W)/W}V_2 $$

cumple que $Z_1,Z_2\stackrel{\text{i.i.d}}{\sim} \text{Normal}(0,1)$.

Cociente de uniformes¶

Del método de Marsaglie sabemos que si $Z_1,Z_2\stackrel{\text{i.i.d}}{\sim} \text{Normal}(0,1)$ entonces podemos expresar

$$ \frac{Z_1}{Z_2}=\frac{\sqrt{-\log(W)/W}V_1}{\sqrt{-\log(W)/W}V_2}=\frac{V_1}{V_2} $$

como un cociente de Uniformes $V_1,V_2$. Podemos generalizar esta idea al considerar, en vez del disco unitario como en el método de Marsaglia, una región

$$ \Omega(f)=\left\{ (v_1,v_2)\, : \, 0\ < v_1 \leq \sqrt{ h(v_2/v_1) } \right\} $$

con $h=cf$ tal que $f$ es la densidad que nos interesa simualr. Si $(V_1,V_2)$ es uniforme en $\Omega$ entonces

$$ X=\frac{V_2}{V_1}\sim \mathcal{L}(f). $$

Podemos nuevamente usar aceptación rechazo si consideramos una rectángulo $R$ tal que $\Omega \subseteq R$. En particular podemos toma $R=(0,a)\times (b_- , b_+)$ con

\begin{align*} a &= \left( \sup\left\{ h(x) \, : \, x\in\mathbb{R} \right\} \right)^{1/2} \\ \\ b_+ &= \left( \sup\left\{ x^2 h(x) \, : \, x\geq 0 \right\} \right)^{1/2} \\ \\ b_- &= -\left( \sup\left\{ x^2 h(x) \, : \, x\leq 0 \right\} \right)^{1/2} \end{align*}

Algoritmo de Ziggurat¶

Drawing

In [192]:
plot(p)
Out[192]:

Drawing

In [240]:
plot(p)
Out[240]:

Denotamos a los rectángulos formados en el ziggurat, empezando de arriba hacia abjo, como $R_1,R_2,\ldots, R_7$ y a la base del ziggurat la denotamos $R_0$, observemos que este último conjunto está conformada por la unión de un rectángulo con la cola de la densidad sin normalizar.

In [205]:
plot(p)
Out[205]:

Algoritmo de Ziggurat¶

1) Elija uniformemente $i\in \left\{ 0,1,\ldots, N-1 \right\}$ y genere $U_1,U_2\sim \text{Uniforme}(0,1)$ independientes. Sea $x=U_1 x_i$

2) Si $$ x < x_{i-1}$$ aceptamos $X=x$ y terminamos; en caso contrario proceder a 3).

3) Si $i=0$ regresamos $X$ simulado de la cola $[x_{N-1},\infty)$ y terminamos; en caso contrario proceder a 4)

4) Si $ U_2 < \frac{f(x)-f(x_i)}{f(x_{i-1})-f(x_i)} $ regrese $X=x$; en caso contrario regrese a 1)

Algoritmo de Ziggurat con operaciones de punto flotante para $N=2^8=256$¶

Para cada $i\in\left\{ 1,2,\ldots, 255\right\}$, defina $k_i = \lfloor 2^{32} x_{i-1}/x_i \rfloor$ y $w_i = x_i/2^{32}$. Para $i=0$ definimos $k_0 = \lfloor x_{255}f(x_{255})/v \rfloor$.

1) Sea $j$ un entero de 32 bits sin signo uniforme y $U\sim \text{Uniforme}(0,1)$. Sea $x= j w_i$

2) Si $$ j < k_i $$ aceptamos $X=x$ y terminamos; en caso contrario proceder a 3).

3) Si $i=0$ regresamos $X$ simulado de la cola $[x_{255},\infty)$ y terminamos; en caso contrario proceder a 4)

4) Si $ U_2 < \frac{f(x)-f(x_i)}{f(x_{i-1})-f(x_i)} $ regrese $X=x$; en caso contrario regrese a 1)

Donde hemos denotado $v$ como el área de los rectángulos $R_0,R_1,\ldots R_{N-1}$. Observemos que

$$ v=x_i(f(x_{i-1})-f(x_i)),\text{ para } i\in\left\{1,2,\ldots , N-1 \right\}, $$

con $x_0=0$, y

$$ v = rf(r) + \int_r^\infty e^{-x^2/2}\big/\sqrt{2\pi}. $$

Esto puede ser calculado numéricamente como se describe a continuación.

Consideramos una función $z(r)$ dada por tomar $v = rf(r) + \int_r^\infty e^{-x^2/2}\big/\sqrt{2\pi}$,

$$ x_i = f^{-1}\left( v/x_{i+1} + f(x_{i+1}) \right) $$

y regresar $v-x_1+x_1f(x_1)$. Usando un método numérico, por ejemplo el método de la bisección, podemos encontrar $z(r)=0$ para fijar el punto de corte $r=x_{N-1}$.

In [235]:
n = 7
function z(r)
    x = r 
    v = r*f(r) + hcubature(f, [r], [100])[1]
    for i in 1:(n-1)
        x = finv(v/x + f(x))
    end
    return v -x + x*f(x)
end
Out[235]:
z (generic function with 1 method)
In [236]:
z(1.1)
Out[236]:
0.9407462535796958
In [237]:
z(3.1)
Out[237]:
-2.1082109732965058
In [238]:
ziggur_r =find_zero(z, (1.1,3.1), Bisection())
Out[238]:
2.338371698247253
In [247]:
R_vec = [ x[i+1]*(f(x[i])-f(x[i+1])) for i in 1:7 ]
Out[247]:
7-element Array{Float64,1}:
 0.17617364011877823
 0.17617364011877745
 0.17617364011877737
 0.17617364011877745
 0.1761736401187774 
 0.17617364011877745
 0.17617364011877754
In [248]:
ziggur_v
Out[248]:
0.17617364011877745
In [372]:
@time BoxMullernormal = BoxMullerNormal(10000);

  0.002312 seconds (10.01 k allocations: 1015.859 KiB)
In [370]:
@time MarsagliaPolarnormal = MarsagliaPolarNormal(10000);

  0.003384 seconds (25.41 k allocations: 2.402 MiB)
In [369]:
@time CocienteUnifnormal = CocienteUnifNormal(10000);

  0.001777 seconds (13.69 k allocations: 1.330 MiB)
In [368]:
@time Zigguratnormal = randn(10000);

  0.000278 seconds (6 allocations: 78.359 KiB)
In [371]:
@time ARSnormal = adaptive_accept_reject([-1.1,1.1],10000,h,dh);

  0.400919 seconds (12.65 M allocations: 304.796 MiB, 13.58% gc time)
In [375]:
plot(p)
Out[375]:
In [503]:
plot(p)
Out[503]:
In [378]:
plot(p)
Out[378]:
In [380]:
plot(p)
Out[380]:
In [383]:
plot(p)
Out[383]:

Generación de variables aleatorias Gamma¶

Algoritmo GD de Ahrens y Dieter (Generating gamma variates by a modified rejection technique, 1982)¶

La idea principal es considerar $X\sim \text{Gamma}(\alpha,1)$ con $\alpha \geq 1$ y usar la transformación $Y=g(X)=\left( \sqrt{\alpha-1/2}+X/2\right)^2$ que hace que $Y$ tenga densidad cercana a la Normal$(0,1)$.

In [505]:
plot(p)
Out[505]:
In [507]:
plot(p)
Out[507]:

Algoritmo MT de Marsaglia y Tsang (Generating gamma variates by a modified rejection technique, 1982)¶

La idea principal es usar aceptación-rechazo para generar $X$ con densidad de probabilidad proporcionala a $e^{g(x)}\leq e^{-x^2/2}$ tal que existe una tranformación $h(X)=Y\sim \text{Gamma}(\alpha,1)$. Este método es usualmente usado con $\alpha\in(1,2)$ junto con el hecho de que para $\hat{\alpha}\in(0,1)$ $\gamma_{\hat{\alpha}+1}\sim\text{Gamma}(\hat{\alpha}+1,1)\implies U^{1/\hat{\alpha}}\gamma_{\hat{\alpha}+1}\sim \text{Gamma}(\hat{\alpha},1)$ para $U\sim \text{Uniforme}(0,1)$. Por lo que el algoritmo puede ser usado para simular variables aleatorias gamma con $\alpha \in(0,2)$.

In [510]:
plot(p)
Out[510]:

Gaussiana Multivariada¶

Un vector Gaussiano multivariado $\pmb{X}=(X_1,\ldots ,X_n)'$ está determinado por su vector de medias $\pmb{\mu}=\mathbb{E}[\pmb{X}]$ y la matriz de varianza-covarianza $\pmb{\Sigma}$ tal que su entrada $(i,j)$ está dada por $\Sigma[i,j]=\text{Cov}\left( X_i, X_j \right)$. Denotamos $\Sigma=\text{Cov}(\pmb{X})$ y más en general $\text{Cov}(\pmb{X},\pmb{Y})$ como la matriz con entradas $\text{Cov}\left( X_i, Y_j \right)$. Observemos que de la linearidad de la covarianza se sigue que para $A$ una matrix $m\times n$ dimensional se tiene que $\text{Cov}(A\pmb{X},\pmb{Y})=A\text{Cov}(\pmb{X},\pmb{Y})$ y $\text{Cov}(\pmb{X},A\pmb{Y})=A'\text{Cov}(\pmb{X},\pmb{Y})$.

Sabemos que si $X_1$ y $X_2$ son variable aleatoriales normales entonces $\text{Cov}(X_1,X_2)=0\implies X_1 \perp X_2$ por lo que es sencillo simular $\pmb{X}$ con entradas Gaussianas tal que $\Sigma=\text{Cov}(\pmb{X})=I$. Por otro lado observamos que si nos intersa simular un vector con entrdas Gaussianas que tenga media $\pmb{\mu}$ y matriz de varianza-covarianza $\Sigma=CC'$ podemos generar $\pmb{Z}$ con entradas Gaussianas tal que $\mathbb{E}[\pmb{Z}]=\pmb{0}$ y $\text{Cov}(\pmb{Z})=I$ con lo que $\pmb{X}=C\pmb{Z}+\pmb{\mu}$ tiene media $\pmb{\mu}$ y

$$ \text{Cov}(\pmb{X})= \text{Cov}(C\pmb{Z})=\text{Cov}(C\pmb{Z},C\pmb{Z})=CC´\text{Cov}(Z)=CC'I=CC'=\Sigma $$
In [141]:
p2
Out[141]:
In [43]:
plot(p1,p2)
Out[43]:
In [145]:
# Consideramos la matriz de Varianza-Covarianza
Σ = [ 1.0  -0.5 ; -0.5 0.75]
Out[145]:
2×2 Array{Float64,2}:
  1.0  -0.5 
 -0.5   0.75
In [147]:
# Con factorización de Cholesky Σ = AA' dada por
A = cholesky(Σ)
Out[147]:
LinearAlgebra.Cholesky{Float64,Array{Float64,2}}
U factor:
2×2 LinearAlgebra.UpperTriangular{Float64,Array{Float64,2}}:
 1.0  -0.5     
  ⋅    0.707107
In [148]:
A.L*A.U
Out[148]:
2×2 Array{Float64,2}:
  1.0  -0.5 
 -0.5   0.75
In [152]:
plot(p)
Out[152]:
In [154]:
p2
Out[154]:
In [56]:
plot(p1,p2)
Out[56]:

Distribuciones Multivariadas Continuas¶

Recordemos que si $X\sim \mathcal{L}(F)$ con $F$ una función cumulativa de probabilidad continua entonces $F(X)\sim \text{Uniforme}(0,1)$. Al considerar vectores de variables aleatorias, esto nos motiva a trabajar con vectores de variables aleatorias marginalmente uniformes. Definimos una cópula $C$ como la función cumulativa de distribución de $\pmb{U}=(U_1,U_2,\ldots , U_d)\in[0,1]^d$ tal que cada $U_i\sim \text{Uniforme}(0,1)$, $i\in\left\{1,\ldots ,d\right\}$. En general, para un vector aleatorio $\pmb{Y}=(Y_1,\ldots , Y_d)$ con funciones cumulativas de probabilidad marginales $F_i$ definimos la cópula de $\pmb{Y}$ como la cópula del vector $(F_1(Y_1),\ldots , F_d(Y_d))$.

Observemos que podemos utilizar la cópula asociada a un vector aleatorio con entradas dependientes para construir vectores con distribuciones marginales que sean distintas pero que mantengan la estructura de dependencia dada por la cópula del vector inicial.

Por ejemplo, si $\pmb{X}\sim \text{Normal}(\pmb{1},)$ entonces $F_i(x)=\phi\left( \frac{x-\mathbb{E}[X_i]}{ \sqrt{\text{Var}(X_i)}}\right)$ con $\phi$ la función cumulativa de probabilidad de una variable aleatoria $Z\sim \text{Normal}(0,1)$. Así las cosas podemos simular $\pmb{U}=(F_1(X_1),\ldots , F_d(X_d))$ y llamamos a la cópula asociada $\textit{cópula Gaussiana}$.

In [156]:
plot(p)
Out[156]:
In [158]:
p2
Out[158]:
In [159]:
p1=histogram(cdf.( Normal(), (A.L*Z)[1,:]/Σ[1,1]^0.5 ), nbins=50, normalize=true, aspect_ratio=1,label="")
title!("Histograma \$U_1\$\n ($(n) observaciones)")
p2=histogram(cdf.( Normal(), (A.L*Z)[2,:]/Σ[2,2]^0.5), nbins=50, normalize=true, aspect_ratio=1,label="")
title!("Histograma \$U_2\$\n ($(n) observaciones)")
plot(p1,p2)
Out[159]:
In [161]:
plot(p)
Out[161]:
In [166]:
p2
Out[166]:
In [168]:
plot(p1,p2)
Out[168]:
In [171]:
plot(p)
Out[171]:
In [173]:
plot(p2)
Out[173]:
In [68]:
plot(p1,p2)
Out[68]:
In [176]:
plot(p)
Out[176]:
In [181]:
plot(p2)
Out[181]:
In [183]:
plot(p1,p2)
Out[183]:
In [185]:
plot(p)
Out[185]:
In [188]:
plot(p2)
Out[188]:
In [196]:
plot(p)
Out[196]:
In [203]:
plot(p2)
Out[203]:
In [205]:
plot(p)
Out[205]:
In [207]:
plot(p2)
Out[207]:
In [209]:
plot(p)
Out[209]:
In [212]:
plot(p2)
Out[212]:
In [214]:
plot(p)
Out[214]:
In [216]:
plot(p2)
Out[216]:
In [218]:
plot(p)
Out[218]:
In [220]:
plot(p2)
Out[220]:
In [222]:
plot(p)
Out[222]:
In [225]:
plot(p2)
Out[225]:

Teorema de Sklar¶

Denotemos $F(x_1,\ldots , x_d)=\mathbb{P}[X_1 \leq x_1, \ldots , X_d \leq x_d]$ la función de distribución cumulativa de $\pmb{X}=(X_1, \ldots , X_d)$, $F_i(x)=\mathbb{P}[X_i\leq x]$ la función cumulativa de probabilidad marginal de $X_i$ para $i\in\left\{1,\ldots , n\right\}$ y $C$ la cópula asociada $\pmb{X}$. Entonces

$$ F(x_1, \ldots , x_d) = C(F_1(x_1), \ldots , F_d(x_d)). $$

Si $F$ tiene densidad $f$ y $F_i$ tiene densidad $f_i$ entonces la densidad de $C$ está dada por

$$ c(u_1,\ldots ,u_d)=\frac{\partial^d}{\partial u_1 \cdots \partial u_d}C(u_1,\ldots , u_d) =\frac{f\left( F_1^{-1}(u_1), \ldots , F_d^{-1}(u_d) \right)}{ f_1\left(F_1^{-1}(u_1)\right)\cdots f_d\left(F_d^{-1}(u_d)\right) }. $$

Observemos que dada la función cumulativa conjunta $F$ y las márginales $F_1,\ldots , F_d$ podemos obtener la cópula asociada $C$ como

$$ C(u_1,\ldots, u_d) = F\left( F_1^{-1}(u_1), \ldots , F_d^{-1}(u_d) \right) $$

Si tenemos un vector aleatorio bivariado $(X_1,X_2)$ con densidad $f_{X_1,X_2}(x_1,x_2)$, marginales $f_{X_1}(x)$, $f_{X_2}(x)$ y condicionales $f_{X_1|X_2=x_2}(x_1)=f_{X_1,X_2}(x_1,x_2)/f_{X_2}(x_2)$, $f_{X_2|X_1=x_1}(x_2)=f_{X_1,X_2}(x_1,x_2)/f_{X_1}(x_1)$. Podemos generar $(X_1,X_2)$ de la siguiente manera

1) Simulamos $X_1\sim \mathcal{L}(f_{X_1})$

2) Simulamos $X_2 \sim \mathcal{L}(f_{X_2|X_1=x_1})$

Si $C_{U_1,U_2}$ es la cópula asociada a $f_{X_1,X_2}$ entonces la función cumulativa de probabilidad asociada a $U_2|U_1$ está dada por

$$\mathbb{P}[U_2\leq u_2|U_1=u_1]=\frac{\partial }{\partial u_1}C(u_1,u_2)=C_{u_1}(u_2)$$

.

Por lo que se puede usar el siguiente algoritmo de simulación.

1) Genere variables aleatorias independientes $V_1,V_2\sim \text{Uniforme}(0,1)$.

2) Tome $U_1=V_1$ y $U_2=C^{-1}_{U_1}(V_2)$.

3) Regrese $X_1=F_1^{-1}(U_1)$ y $X_2=F_2^{-1}(U_2)$

Ejemplos de Cópulas¶

1) Cópula de Clayton. Sea $\theta> 0$

$$ C(u_1,\ldots,u_d)=\left( \sum_{i=1}^d u_i^{-\theta}-d+1 \right)^{-\frac{1}{\theta}}. $$

2) Cópula de Gumbel. Sea $\theta > 1$

$$ C(u_1,\ldots,u_d)=\exp\left( - \left( \sum_{i=1}^d -\log(u_i)^{\theta} \right)^{\frac{1}{\theta}} \right). $$

3) Cópula de Frank. Sea $\theta > 0$

$$ C(u_1,\ldots,u_d)=-\frac{1}{\theta}\log\left( 1+ \frac{\prod_{i=1}^d\left(e^{-\theta u_i}-1\right)}{\left(e^{-\theta}-1\right)^{d-1}} \right). $$
In [286]:
p
Out[286]:
In [290]:
p2
Out[290]:
In [293]:
p
Out[293]:
In [296]:
p2
Out[296]:
In [299]:
plot(p)
Out[299]:
In [301]:
p2
Out[301]: